Cool things to do in Mathematica

Mit Mathematica kann man viele coole Sachen anstellen. Ich hab hier mal ein paar zusammengefasst um euch einen Eindruck zu vermitteln wozu dieses Programm alles fähig ist.
Alle Beispiele sind ein bisschen kompliziert und mathematisch/programmiertechnisch manchmal recht anspruchsvoll aber ich habe mich bemüht,
alles so gut wie möglich zu erklären/kommentieren.

“Vielleicht ist nächstes Jahr auch jemand von euch im Tutorium und fühlt sich animiert, ebenfalls eine Mathematica Einführung zu machen und dieses Notebook evtl. um weitere coole Sachen zu ergänzen.”
  -Julien Kluge 2015

Fourierzerlegung der Wand des Gerthsen-Hörsaals

Hast du dich schonmal gefragt, wieso die Wand des Gerthsen-Hörsaals (Physik Institut - Newtonstraße 15.) ein so schreckliches Muster ziert? Dieses scheinbar missglückte Experiment von ‘Kunst am Bau’ hat aber tatsächlich Hintergründe. 
Wer jetzt nicht gespoilert werden möchte ließt nicht weiter...
Es handelt sich nämlich um die Fouriertransformation des HU-Logos. Wem die Fourierzerlegung nichts sagt, der kann mal im Internet nach Literatur suchen. Für alle anderen:
die Fourierzerlegung kann nicht nur über Gleichungen sondern auch diskret über Daten gemacht werden. Ich zeige hier ein Beispiel wie wir die Wand vor dem Gerthsen entschlüsseln können.

gerthsenWallImage=; (*Laden des Bildes der Wand in eine Variable*)

gerthsenWallData=ImageData[ColorConvert[gerthsenWallImage,"Grayscale"]][[All,All,1]]; (*Transformieren des Bildes zu Schwarz-Weiß und anschliessende Transformation in eine Datenmatrix*)

gerthsenWallFourierData=Abs[Fourier[gerthsenWallData]]^2;
(*Durchführung der FFT (Fast Fourier Transformation) und Bildung des Quadrates des Absolutbetrages, da die Daten komplex sein können und das Quadrat um die Schärfe zu erhöhen*)

gerthsenWallFourierImage=Image[gerthsenWallFourierData]
(*Wiederherstellung des Bildes aus den Fourier transformierten Daten*)

gerthsenWallFourierImageDim=ImageDimensions[gerthsenWallFourierImage]; (*speichern der Bildgröße*)
ImageTrim[ImageCrop[gerthsenWallFourierImage,gerthsenWallFourierImageDim*2,Padding->"Periodic"],{{0,0},gerthsenWallFourierImageDim}]
(*Vergrößern und anschliessendes ausschneiden des Bildes um es schön zu zentrieren*)

Viele Wege führen nach Rom(/zum Ergebnis)

Mathematica ist eine multiparadigmen-Sprache. Das bedeutet, sie ist in der Lage, auf viele verschiedene Art und Weisen Aufgaben zu erledigen. Das ist einerseits gut, andererseits schlecht. Gut, weil man immer die effizienteste Art und Weise wählen kann etwas zu tun oder die am kürzesten zu schreibende. Schlecht, weil Redundanz auch sehr leicht zu Problemen führen kann, es manchmal schmerzhaft ist, nachzutesten was denn nun die schnellste Form der Implementierung ist und es den Code auch recht schnell sehr unübersichtlich machen kann.
(Wenn man übrigens nach einem “Hauptparadigma” sucht, welches Mathematica definiert, dann ist dieses wohl die Regelbasiertheit auch wenn man davon so gut wie nie was mitbekommt...)
Im folgenden führe ich mal einige Wege auf wie man Zahlen von eins bis n aufsummieren kann.

(1) Die einfachste Möglichkeit ist es wohl, die Zahlen einfach einzugeben...wird nur blöd für große n

1+2+3+4+5+6+7+8+9+10

(2, 3) Man denke auch daran, dass man einfach auf english fragen kann:

(4, 5, 6) Wenn ihr schon Erfahrung mit anderen Programmiersprachen habt, dann ist der Prozedurale Ansatz wohl das, womit ihr am meisten vertraut seid. (Mathematica ist hier übrigens relativ langsam)

var1=0;Do[var1+=i,{i,1,10}];var1 (*Procedural*)

var2=0;For[i=1,i≤10,++i,var2+=i];var2 (*Procedural*)

var3=0;var4=0;While[var4≤9,var3+=(var4+=1)];var3

(7) Oder dann doch lieber über Rekursion?

g[n_]:=If[n==1,1,g[n-1]+n];g[10]

(8) Jeder richtige Mathematiker unter euch würde mir natürlich sagen, dass die Summe von 1 bis N eine Identität hat:

(10+1)*10/2

(9) Wir können aber natürlich auch sofort die Sum-Funktion anstrengen:

Sum[i,{i,1,10}] (*Declarative*)

(10) Und das lässt sich auch über Palettes direkt eingeben:

(∑_(i=1)^10)i

(11, 12) Wir können die Identität natürlich auch über Funktionen callen oder über Regeln einsetzen

f[n_]:=(n+1)*n/2;f[10] (*Functional*)

((n+1)*n/2)/.n->10 (*Symbolic, Rules Based*)

(13, 14, 15) Es gibt wirklich viele Ansätze das Listenbasiert zu lösen über: Range, Table, Total & Sum

list=Range[10];Sum[list[[i]],{i,1,Length[list]}] (*Procedural, Declarative, List Based*)

Total[Table[i,{i,1,10}]] (*Declarative, List Based*)

Total[Range[10]] (*Declarative, List Based*)

(16, 17, 18, 19 , 20) Und dem entsprechend auch über Funktionaloperatoren mit Range, Table & Total

Apply[Plus,Table[i,{i,1,10}]] (*Declarative, Functional Operators, List Based*)

Apply[Plus,Range[10]] (*Declarative, Functional Operators, List Based*)

Plus@@Table[i,{i,1,10}] (*Declarative, Functional Operators, List Based*)

Plus@@Range@10 (*Declarative, Functional Operators, List Based*)

Range@10//Total (*Declarative, Functional Operators, List Based*)

(21, 22) Man kann es auch mit den Regeln übertreiben

addr[10]//.{addr[1]->1,addr[n_Integer]:>n+addr[n-1]} (*Rules Based*)

Merge[Table[a->i,{i,1,10}],Total][[1]] (*Functional Operators, List Based, Rules Based*)

(23) Oder generell übertreiben...

list=Range[10];While[Length[list]>1&&list[[-1]]≠0,If[Mod[Length[list],2]≠0,list=list~Join~{0}];list=Partition[list,2]/.{x_,y_}->(x+y)];list[[1]] (*Procedural, List Based, Rules Based, Pattern Matching*)

(24, 25) Man kann Mathematica natürlich auch zu einer Scriptsprache reduzieren

ToExpression["1+2+3+4+5+6+7+8+9+10"]

ToExpression[StringReplace[ToString[Table[i,{i,10}]],{","->"+"," "->"","{"->"","}"->""}]] (*Interpreted, List Based*)

(26, 27 , 28, 29) Oder man fängt an zu cheaten und summiert nicht mehr auf sondern sucht Funktionen die das ganze Beschreiben können und wertet diese dann aus:

Round[(LinearModelFit[{1,3},{x,x^2},x,IncludeConstantBasis->False]["BestFit"])/.x->10] (*Rules Based, OOP*)

Round[(NonlinearModelFit[{1,3},a*x+b*x^2,{a,b},x]["BestFit"])/.x->10] (*Rules Based, OOP*)

Round[Interpolation[{1,3,6,10},Method->"Spline"][10]] (*Declarative, Functional*)

FindSequenceFunction[{1,3,6,10}][10] (*Declarative, Functional*)

(30) Oder man versucht “künstliche Lernalgorithmen” (Ich hütte mich vor dem Ausdruck: ‘künstliche Intelligenz’) zu benutzen (auch wenn die nicht wirklich genau ist...)

Predict[Flatten[{#->Total[Range[#]]}&/@{1,2,3,4,5,6,7,8,9,11}],Method->"NearestNeighbors"][10]

30 Möglichkeiten für ein und die Selbe Sache ... nicht schlecht.
To be continued ...

Machinelles Lernen

Mathematica stellt für machinelles lernen zwei Funktionen zur Verfügung ‘Predict’, ‘Classify’. Beide müssen mit einem Set trainiert werden und können nach dem Training befragt werden.
In diesem Fall stell ich euch mal zwei Möglichkeiten vor.

Hättest du das Titanic - Unglück überlebt?

Als erstes laden wir mal die Überlebensdaten/Sterbedaten von 1309 Passagieren der Titanic von Wolfram-Research herunter.

titanicData=ExampleData[{"MachineLearning","Titanic"},"Data"];

Die Daten kommen in der Form {Klasse,Alter,Geschlecht}->(überlebt/gestorben)
{{“1st”, “2nd”, “3rd”}, {0.0-100.0, Missing[]}, {“female”,”male”}}->{survived, died}
Schauen wir uns mal 5 Beispiele an:

RandomSample[titanicData,5]//TableForm

titanicClassify=Classify[titanicData,Method->"LogisticRegression"] (*trainieren wir unseren Algorithmus*)

Zum Zeipunkt an dem ich das hier schreibe, bin ich : männlich, 20.786 Jahre alt und hab soviel Geld das es wohl nur für ein Ticket dritter Klasser reicht ..
Mal sehen wie es für mich ausgesehen hätte, wäre ich jetzt auf der Titanic:

titanicClassify[{"3rd",20.789,"male"}]

fuck ... 

Naja schauen wir uns mal die Wahrscheinlichkeiten an, die ich in allen Klassen gehabt hätte:

titanicClassify[{"3rd",20.789,"male"},{"Probability","survived"}]*100"%"
titanicClassify[{"2nd",20.789,"male"},{"Probability","survived"}]*100"%"
titanicClassify[{"1st",20.789,"male"},{"Probability","survived"}]*100"%"

na ganz klasse...
Überlebst du Wahrscheinlicher?

Welchem Autor wärst du am nächsten?

Classify kann ohne Probleme dazu eingesetzt werden, zu bestimmen welche Texte zu welchen Autor gehören.
Als erstes laden wir mal texte von ihnen herunter.

shakespearHamlet=Import["http://www.gutenberg.org/cache/epub/2265/pg2265.txt"];
(*Hamlet von William Shakespear*)
dickensOliverTwist=Import["http://www.gutenberg.org/files/47529/47529-0.txt"];
(*Oliver Twist vol. 1 von Charles Dickens*)
twainMysteriousStranger=Import["http://www.gutenberg.org/files/50109/50109-0.txt"];
(*The Mysterious Stranger von Mark Twain*)
bibleKingJames=Import["http://www.gutenberg.org/cache/epub/10/pg10.txt"];
(*King James Version der Bibel (altes Testament)*)

authorClassify=Classify[{shakespearHamlet->"William Shakespear",dickensOliverTwist->"Charles Dickens",twainMysteriousStranger->"Mark Twain",bibleKingJames->"King James Bibel"}]

Jetzt füge ich da einfach mal ein Text ein, den ich so für den Moodle-Kurs des Brückenkurses 2015 geschrieben hatte:

authorClassify["If you need a competent person for questions regarding, studying abroad in germany, studying physics on the Humboldt University of Berlin, informations for your courses, timetables,..., etc., please feel free to contact us (You can use german, english or swedish). We'll be glad to help you. You can nevertheless contact any other person you wish as long as you are writing in german or english. sincerely"]

Ich denke, ich kann damit leben, so ähnlich wie Mark Twain zu schreiben.

authorClassify["you all shall pass the study of the holy physics"]

Hehe, war ja klar : D

Finanzverläufe vorhersagen (oder auch nicht)

Wir versuchen mal, mithilfe von Regression den Kursverlauf einer Aktie für kleine Zeiten vorherzusagen.

appleStock=FinancialData["NASDAQ:AAPL",{"Jan. 1, 2014",Today}];

DateListPlot[appleStock,PlotRange->{{"Jan. 1, 2014","Jan. 1, 2016"},Automatic}]
(*naja nehmen wir als Versuchskaninchen mal die Apple-Aktie*)

(*Wir müssen die Daten in Zahlen überführen um darüber regressieren zu können*)
appleDayStock=Table[{N[QuantityMagnitude[DateDifference[appleStock[[Length[appleStock]-i+1,1]],Today]]/100,5],appleStock[[i,2]]},{i,Length[appleStock]}];

(*Inwieweit können wir das ganze über Sinus Funktionen annähern um die nächsten 5 Tage vorherzusagen
Achtung: das braucht eine lange Zeit zum rechnen*)
appleLinSinFit=LinearModelFit[appleDayStock[[1;;(Length[appleDayStock]-4)]],Table[Sin[a*x],{a,0.2,3,0.2}]~Join~{1}~Join~Table[Sin[b*x],{b,4,20}],x];
Show[ListPlot[appleDayStock[[1;;(Length[appleDayStock]-4)]],PlotRange->{{0,7.1},{65,140}},PlotTheme->"Detailed",AspectRatio->1/3,PlotMarkers->"x"],ListPlot[appleDayStock[[(Length[appleDayStock]-4);;Length[appleDayStock]]],PlotStyle->Red,PlotMarkers->"O"],Plot[appleLinSinFit["BestFit"],{x,0,7.1},PlotStyle->Green],ImageSize->1100]

ToString[R^2,TraditionalForm]<>"="<>ToString[appleLinSinFit["AdjustedRSquared"]]
ToString[χ^2,TraditionalForm]<>"="<>ToString[appleLinSinFit["ANOVATableMeanSquares"][[-1]]]

Naja immerhin 3 Tage sind in unserem Vorhersageintervall...reich werden wir dadurch nicht
Dieses Modell lässt sich aber weiter sukzessive Verbessern durch Einsatz neuer Funktionen, Gewichtungen und statistischen Vorhersagen.
Hier eine Anfängervariante mit plain Mathematica-Befehlen.

tsm=TimeSeriesModelFit[appleStock,"MA"];
tsf=TimeSeriesForecast[tsm,{50}];
DateListPlot[{appleStock,tsf+3},ImageSize->800,PlotRange->{{"Jan. 1. 2015",All},{100,135}},AspectRatio->1/3]

Nie wiederkehrender Plot

WIr benutzen die RiemannSiegelZ Funktion um einen Plot zu erzeugen, der sich fortlaufend ändert aber niemals der Selbe ist.

c=20;
xMax=5;xMin=-5;
yMax=5;yMin=-5;
Dynamic[c=1.3^(1/c)*c;
dataSiegel=Table[{RiemannSiegelZ[j-c],RiemannSiegelZ[j]},{j,1,50,0.1}];
xMax=Max[dataSiegel[[All,1]],xMax];xMin=Min[dataSiegel[[All,1]],xMin];
yMax=Max[dataSiegel[[All,2]],yMax];yMin=Min[dataSiegel[[All,2]],yMin];
ListPlot[dataSiegel,PlotRange->{{xMin,xMax},{yMin,yMax}},ImageSize->Large]]

Elektrische Felder

Das elektrische Feld mehrerer Ladungen (Ort (OverVector[(Subscript[s, i])])) ist mit Superposition an einem Ort (OverVector[r]) gegeben mit

(OverVector[E])((OverVector[r]))=(1/(4π (Subscript[ϵ, 0])))(∑_(i=1)^N)(((Subscript[q, i]) ((OverVector[r])-(OverVector[(Subscript[s, i])])) )/(|(OverVector[r])-(OverVector[(Subscript[s, i])])((|)^3)))

Wie Professor Busch jetzt sagen würde dieser Vorfaktor, die Physikpolizei, ist erstmal nicht nötig. Die Formel wird also zu:

(OverVector[E])((OverVector[r]))=(∑_(i=1)^N)(((Subscript[q, i]) )/(|(OverVector[r])-(OverVector[(Subscript[s, i])])((|)^2)))(OverHat[(Subscript[s, i])])

Um diese Formel anzuwenden erzeugen wir erstmal zufällige Ladungen.

chargeCount=10; (*Anzahl der Ladungen*)
charges=RandomReal[{1,4},chargeCount]; (*Zufällige Ladung zwischen 1 und 4 Coulomb*)
chargeR=RandomReal[{-10,10},{chargeCount,2}];(*Positionen der Ladungen zwischen +-10 im 2D Raum*)

Der elektrische Feldvektor an einem Punkt {x,y} ist also gegeben mit:

eFeldVec[x_,y_,count_,chargesP_,chargesRP_]:=Sum[chargesP[[i]]/(Norm[{x-chargesRP[[i,1]],y-chargesRP[[i,2]]}]^2)*Normalize[{x-chargesRP[[i,1]],y-chargesRP[[i,2]]}],{i,1,count}]

Am Punkt {0,0} ist der Feldvektor also:

eFeldVec[0,0,chargeCount,charges,chargeR]

Wir können nun also es zu plotten (benutzt niemals VectorPlot...):

StreamPlot[eFeldVec[x,y,chargeCount,charges,chargeR],{x,-10,10},{y,-10,10},ImageSize->600]

Zeichnen wir auch mal die Ladungen mit ein

Show[StreamPlot[eFeldVec[x,y,chargeCount,charges,chargeR],{x,-10,10},{y,-10,10}],ListPlot[chargeR,PlotStyle->Red]]

Es ergiebt auch sehr interessante Grenzverläufe der Feldlinien bei regelmäßigen Mustern

chargeCountB=13; (*Anzahl der Ladungen*)
chargesB=Table[1,chargeCountB]; (*Zufällige Ladung zwischen 1 und 4 Coulomb*)
chargeRB=Table[7*{Cos[2Pi/chargeCountB*i],Sin[2Pi/chargeCountB*i]},{i,1,chargeCountB}];(*Positionen der Ladungen zwischen +-10 im 2D Raum*)

Show[StreamPlot[eFeldVec[x,y,chargeCountB,chargesB,chargeRB],{x,-10,10},{y,-10,10}],ListPlot[chargeRB,PlotStyle->Red]]

Besonders dann, wenn man die Symmetrie stört

chargeCountC=11; (*Anzahl der Ladungen*)
chargesC=Table[1,chargeCountC]; (*Zufällige Ladung zwischen 1 und 4 Coulomb*)
chargeRC=Table[7*{Cos[2Pi/(chargeCountC+2)*i],Sin[2Pi/(chargeCountC+2)*i]},{i,1,chargeCountC}];
Show[StreamPlot[eFeldVec[x,y,chargeCountC,chargesC,chargeRC],{x,-10,10},{y,-10,10}],ListPlot[chargeRC,PlotStyle->Red]]

Mit unterschiedlichen Vorzeichen für die Ladungen:

chargeCountD=10; (*Anzahl der Ladungen*)
chargesD=RandomChoice[{-4,-3,-2,-1,1,2,3,4},chargeCountD];
chargeRD=RandomReal[{-10,10},{chargeCountD,2}];
Show[StreamPlot[eFeldVec[x,y,chargeCountD,chargesD,chargeRD],{x,-10,10},{y,-10,10}],ListPlot[chargeRD,PlotStyle->Red],ImageSize->600]

Nun probieren wir das mal im 3 dimensionalen:

chargeCount3D=3; (*Anzahl der Ladungen*)
charges3D=RandomReal[{1,4},chargeCount];
chargeR3D=RandomReal[{-10,10},{chargeCount,3}];(*Positionen der Ladungen zwischen +-10 im 3D Raum*)

eFeldVec3D[x_,y_,z_]:=Sum[charges3D[[i]]/(Norm[{x-chargeR3D[[i,1]],y-chargeR3D[[i,2]],z-chargeR3D[[i,3]]}]^2)*Normalize[{x-chargeR3D[[i,1]],y-chargeR3D[[i,2]],z-chargeR3D[[i,3]]}],{i,1,chargeCount3D}]

VectorPlot3D[eFeldVec3D[x,y,z],{x,-10,10},{y,-10,10},{z,-10,10}]

Und hier schlägt das inverse-quadrat in voller Hässlichkeit zu (weswegen ich auch geschrieben hab: benutz nie VectorPlot [und weil der Plot mit Polstellen nicht umgehen kann]). Aber hier sind wir gezwungen den VectorPlot zu benutzen, da es leider keine Version von StreamPlot in 3D gibt.
Um eine schöne Visualisierung zu haben, lasst und mal kurz die Physik über Bord werfen und machen es linear-invers abhängig.

eFeldVec3DW[x_,y_,z_]:=Sum[charges3D[[i]]/(Norm[{x-chargeR3D[[i,1]],y-chargeR3D[[i,2]],z-chargeR3D[[i,3]]}])*Normalize[{x-chargeR3D[[i,1]],y-chargeR3D[[i,2]],z-chargeR3D[[i,3]]}],{i,1,chargeCount3D}]

VectorPlot3D[eFeldVec3DW[x,y,z],{x,-10,10},{y,-10,10},{z,-10,10},VectorScale->Medium]

Viel besser.

Berechnung des Mehrteilchenproblems

Die Berechnung mehrer Teilchen mit nicht vernachlässigbarer Masse die mit Kräften gegenseitig aufeinander wirken, ist nur für zwei Teilchen analytisch lösbar. Für ein System mit mehr als zwei Teilchen wird mit Differentialgleichungen beschrieben, die nicht mehr (bis jetzt) lösbar sind. Man kann das System dennoch Stückweise berechnen. Das tun wir hier.

pCount=5; (*Anzahl an Teilchen*)
dt=0.00001; (*Teilstück der Zeit*)
pPos=Table[{Sin[2Pi/pCount*i],Cos[2Pi/pCount*i]},{i,pCount}]; (*Positionen der Teilchen*)
pVel=0.01/0.00001*Table[{Sin[2Pi/pCount*i+Pi/2],Cos[2Pi/pCount*i+Pi/2]},{i,pCount}]; (*Initialgeschwindigkeit der Teilchen*)
pProp=RandomInteger[{1,5},pCount]; (*Zufällig gewählte Größe der anziehenden Eigenschaft (z.B. Masse, Ladung, etc. ...)*)
camPos={0,0}; (*Position der Kamera*)
Dynamic[
pVel=Table[pVel[[i]]+
Sum[If[i==j,{0,0},Normalize[pPos[[j]]-pPos[[i]]]*(Product[pProp[[k]],{k,pCount}])/Max[0.1,(Norm[pPos[[i]]-pPos[[j]]]^2)]],{j,pCount}]/pProp[[i]]
,{i,pCount}]; (*Berechnung der neuen Geschwindigkeiten der Teilchen*)
pPos+=dt*pVel; (*Berechnung des neuen Ortes über die Zeit t*)
camPos={Mean[pPos[[All,1]]],Mean[pPos[[All,2]]]}; (*Berechnung der neuen Kameraposition*)
ListPlot[pPos,ImageSize->600,PlotRange->{{camPos[[1]]-5,camPos[[1]]+5},{camPos[[2]]-5,camPos[[2]]+5}},PlotTheme->"Scientific"] (*plotting*)
]

Das lässz sich auch sehr gut mit dem letzten Beispiel (eketrische Feldlinien plotten) verbinden:

qCount=10;
qPositions=RandomInteger[{-8,8},{qCount,2}];
qCharges=RandomChoice[{-4,-3,-2,-1,1,2,3,4},qCount];
eleP={0,0};
eleV={0,0};(*5*RandomInteger[{2,5},2];*)
eleQ=2;
baseBackGraph=Show[ListPlot[qPositions,PlotRange->{{-10,10},{-10,10}},PlotTheme->"Scientific",AspectRatio->1],StreamPlot[eFeldVec[x,y,qCount,qCharges,qPositions],{x,-10,10},{y,-10,10}]];
Dynamic[
eleV=eleV+Sum[Normalize[eleP-qPositions[[i]]]*eleQ*qCharges[[i]]/(Max[0.3,Norm[eleP-qPositions[[i]]]])^2,{i,qCount}];
eleP+=0.005*eleV;
Show[baseBackGraph,ListPlot[{eleP},PlotRange->{{-10,10},{-10,10}},PlotStyle->Red],ImageSize->600]
]

Auch wenn die entstehenden DGL’s nicht gelöst werden können, sind wir in der Lage sie numerisch zu berechnen. Mathematica hat diese Fähigkeiten eingebaut mit NDSolve.

qCount2=15;
qPositions2=RandomReal[{-1,1},{qCount2,2}];
qCharge2=RandomChoice[{-1,1,2,3,4},qCount2];
tStart2=0;
tEnd2=5;
startVector2={0,0};
startVelocity2=RandomReal[{-2,2},2];
s2={x[t],y[t]}/.First[Quiet[
NDSolve[{
{x''[t],y''[t]}==Sum[qCharge2[[i]]/((Max[0.2,Norm[{x[t],y[t]}-qPositions2[[i]]]])^2)*Normalize[qPositions2[[i]]-{x[t],y[t]}],{i,qCount2}],
{x[0],y[0]}==startVector2,
{x'[0],y'[0]}==startVelocity2
},{x[t],y[t]},{t,tStart2,tEnd2}]
]];
Animate[Show[ParametricPlot[s2,{t,tStart2,u},PlotRange->{{-2,2},{-2,2}}],ListPlot[qPositions2,PlotStyle->Red]],{u,tStart+0.01,tEnd,0.01},AnimationRate->1/2,AnimationRunning->False]

Und das ganze in 3D

qCount3=15;
qPositions3=RandomReal[{-1,1},{qCount3,3}];
qCharge3=RandomChoice[{-1,1,2,3,4},qCount3];
tStart3=0;
tEnd3=5;
startVector3={0,0,0};
startVelocity3=RandomReal[{-2,2},3];
s3={x[t],y[t],z[t]}/.First[Quiet[
NDSolve[{
{x''[t],y''[t],z''[t]}==Sum[qCharge3[[i]]/((Max[0.2,Norm[{x[t],y[t],z[t]}-qPositions3[[i]]]])^2)*Normalize[qPositions3[[i]]-{x[t],y[t],z[t]}],{i,qCount3}],
{x[0],y[0],z[0]}==startVector3,
{x'[0],y'[0],z'[0]}==startVelocity3
},{x[t],y[t],z[t]},{t,tStart3,tEnd3}]
]];
Show[ParametricPlot3D[s3,{t,tStart3,tEnd3},PlotRange->{{-1.2,1.2},{-1.2,1.2},{-1.2,1.2}}],ListPointPlot3D[qPositions3,PlotStyle->Red]]

Wie in diesem Kapitel schon mit der finite Elemente Methode gezeigt, lässt sich das Mehrkörperproblem mit nicht starren Körpern auch über NDSolve lösen:

tEnd=4.5;
sol=solution=NDSolve[{
{x1[0],y1[0]}==0.5*{Cos[2Pi/3],Sin[2Pi/3]},
{x1'[0],y1'[0]}=={0,0.5},
{x2[0],y2[0]}==0.5*{Cos[4Pi/3],Sin[4Pi/3]},
{x2'[0],y2'[0]}=={-0.1,-0.5},
{x3[0],y3[0]}==0.5*{Cos[6Pi/3],Sin[6Pi/3]},
{x3'[0],y3'[0]}=={0.5,-0.5},
{x1''[t],y1''[t]}==1/(Norm[{x2[t]-x1[t],y2[t]-y1[t]}])^2*Normalize[{x2[t]-x1[t],y2[t]-y1[t]}]+1/(Norm[{x3[t]-x1[t],y3[t]-y1[t]}])^2*Normalize[{x3[t]-x1[t],y3[t]-y1[t]}],
{x2''[t],y2''[t]}==1/(Norm[{x1[t]-x2[t],y1[t]-y2[t]}])^2*Normalize[{x1[t]-x2[t],y1[t]-y2[t]}]+1/(Norm[{x3[t]-x2[t],y3[t]-y2[t]}])^2*Normalize[{x3[t]-x2[t],y3[t]-y2[t]}],
{x3''[t],y3''[t]}==1/(Norm[{x1[t]-x3[t],y1[t]-y3[t]}])^2*Normalize[{x1[t]-x3[t],y1[t]-y3[t]}]+1/(Norm[{x2[t]-x3[t],y2[t]-y3[t]}])^2*Normalize[{x2[t]-x3[t],y2[t]-y3[t]}]
},{x1[t],y1[t],x2[t],y2[t],x3[t],y3[t]},{t,0,tEnd}];
s1={x1[t],y1[t]}/.sol;
s2={x2[t],y2[t]}/.sol;
s3={x3[t],y3[t]}/.sol;
Animate[Show[ParametricPlot[s1,{t,0,u},AspectRatio->1,PlotRange->{{-2,2},{-2,2}},PlotStyle->RGBColor[1,0.7,0]],
ParametricPlot[s2,{t,0,u},AspectRatio->1,PlotStyle->RGBColor[0,0.7,1]],ParametricPlot[s3,{t,0,u},AspectRatio->1,PlotStyle->RGBColor[1,0,1]],ListPlot[{s1/.t->u,s2/.t->u,s3/.t->u},PlotStyle->Red],ImageSize->600],{u,0.1,tEnd},AnimationRate->0.1,AnimationRunning->False]

Selbstplottende Formel

Die sogenannte Tupper’s Sel-Referentail Formula, ist eine Ungleichung, die in der Lage ist, jedes Bild welches 17 Pixel hoch ist, als Regionplot darzustellen.
Die Formel hat folgende Form:

(1/2)<⌊mod(⌊(y/17)⌋(2^(-17⌊x⌋-mod(⌊y⌋,17))),2)⌋

Wobei ⌊⌋ die Floor Funktion und mod(m,n) der nte-Modulo ist.
Einen Plot erzeugen wir mit folgendem Code:

tupper1=Table[If[0.5<Floor[Mod[Floor[(y)/17]*2^((-17)*Floor[x]-Mod[Floor[(y)],17]),2]],1,0],{y,0,16},{x,0,100}];
MatrixPlot[tupper1]

Noch nicht viel zu sehen. (Das macht aber auch Sinn wie wir gleich herrausfinden werden.)
Wenn wir der y-Komponente, nun einen (kleines) Offset (k) geben aber, erhalten wir folgenden Plot:

k=;
(*klein ist relativ...*)
tupper2=Table[If[0.5<Floor[Mod[Floor[(y+k)/17]*2^((-17)*Floor[106-x]-Mod[Floor[(y+k)],17]),2]],1,0],{y,0,16},{x,0,106}];
MatrixPlot[tupper2,ImageSize->Large]

Und hier die Überraschung. Diese riesige Zahl k als Offset für die y-Komponente, führt dazu das sie Formel sich selbst plottet!
Es stellt sich herraus, dass sich für ein beliebiges 1bit-Bitmap Bild die konstante Berechnen lässt, mit der die Formel  dieses Bild plotten lässt.
Der Prozess zur Berechnung dieser Konstante ist erschreckend einfach. Man weise jedem Bildpunkt der eingefärbt werden soll eine eins zu, allen anderen eine Null. Nun reiht man diese Zahlen Spaltenweise hintereinander, multipliziert diese Zahl (als Binär gesehen) mit der Zahl 17 und rechnet das dann noch in seine Dezimaldarstellung um. Fertig.
Wir versuchen uns mal selbst daran.

tupperImg=;

Wir wandeln das Bild in die korrespondierenden Daten um.

tupperImgData=ImageData[Binarize[tupperImg],"Bit"];

Jetzt folgt der Prozess zur Berechnung der konstante.

cK=ToExpression["2^^"<>StringJoin[(StringJoin/@Map[ToString,Transpose[tupperImgData],{2}])]]*17

Auf gehts plotten:

tupper3=Table[If[0.5<Floor[Mod[Floor[(16-y+cK)/17]*2^((-17)*Floor[498-x]-Mod[Floor[(16-y+cK)],17]),2]],1,0],{y,0,16},{x,0,498}];
MatrixPlot[tupper3,ImageSize->Large]

Wie die vielen 17nen erahnen lassen, lässt sich das erweitern wenn man möchte:

huLogo=;
ImageDimensions[huLogo]

Es folgt:

(1/2)<⌊mod(⌊(y/100)⌋(2^(-100⌊x⌋-mod(⌊y⌋,100))),2)⌋

Die Berechnung wird genauso ausgeführt wie zuvor nur mit 100 statt 17 multipliziert:

huLogoTupperData=ImageData[Binarize[huLogo],"Bit"];
hK=ToExpression["2^^"<>StringJoin[(StringJoin/@Map[ToString,Transpose[huLogoTupperData],{2}])]]*100

Diese huebsche Zahl stellt das Offset da, mitdem wir das HU-Logo als Offset in der Funktion finden.
also: plotten!

tupper4=Table[If[0.5<Floor[Mod[Floor[(y+hK)/100]*2^((-100)*Floor[x]-Mod[Floor[(y+hK)],100]),2]],1,0],{y,100,1,-1},{x,100,1,-1}];
MatrixPlot[tupper4,ImageSize->Large]

Wer jetzt übrigens glaubt, dass das ein super Weg wäre binäre Bitmaps zu komprimieren schneidet sich gewaltig. Braucht das Original-Bitmap genau 10000 Bits um komplett unkomprimiert gespeichert zu sein, benötigt die Zahl, die unser Offset darstellt, genau so viele Bits:

StringLength[ToString[BaseForm[hK,2]]]

und ist damit gute 48% größer.

Eigenfunktionen des LaPlace-Operators

Der LaPlace-Operator ∆ ist das quasi-Quadrat des Nabla-Operators ∇. 
Für diesen lassen sich Funktionen ϕ finden, die folgende Gleichung erfüllen:

-∆ϕ = λ·ϕ

für ein beliebiges λ∈ℝ\{0} . Die sogenannten Eigenfunktionen des Operators.
Diese stellen sehr hübsche Funktionen im Raum da, die man mal gesehen haben sollte.

(Um Rechenzeit zu sparen, wird das ganze Multithreaded berechnet - ParallelTable)

Ω=Ball[{0,0,0},2]; (*Region einer Kugel mit Radius 2*)
sampleCount=20; (*wie viele Funktionen sollen berechnet werden. Immer: ein vielfaches von 5*)
colorDataListDPlot={Blue,Green,Red,Yellow};
{ϕValues,ϕFunctions}=DEigensystem[{-Laplacian[ϕ[x,y,z],{x,y,z}],DirichletCondition[ϕ[x,y,z]==0,True]},ϕ[x,y,z],{x,y,z}∈Ω,sampleCount];
Grid[Partition[ParallelTable[DensityPlot3D[Evaluate[N[ϕFunctions[[i]]]],{x,y,z}∈Ω,ColorFunction->(Blend[colorDataListDPlot,#]&),ImageSize->200],{i,sampleCount}],Min[sampleCount,5]],Frame->All]